#include "bc_orthogonalDischarge_velocity.hpp"
#include "geometry.hpp"
#include "exceptions.hpp"
#include "cubicRoot.hpp"
#include "field.hpp"
#include "sweFunctions.hpp"


namespace hms {

using std::string_literals::operator""s;

/* positive qBoundary => inflow  */
/* negative qBoundary => outflow */

#define CHEAD(RETURN_TYPE) \
template<typename MeshType> \
RETURN_TYPE OrthogonalDischarge_velocity<MeshType>::


CHEAD() OrthogonalDischarge_velocity(
	Field<MeshType>& v,
	const BoundaryPatch& p,
	Field<MeshType>& dFieldArg,
	scalar qBoundaryArg
) :
	BaseType(Type::orthogonalDischarge_velocity, v, p),
	m_transformer {v, p}
{
	qBoundary(qBoundaryArg);
	dField( dFieldArg );
	m_calcBuf.resize( p.size(), bufCols );
}

CHEAD(auto) qBoundary() const -> scalar {
	return m_qBoundary;
}

CHEAD(void) qBoundary(scalar qBoundaryArg){
	m_qBoundary = qBoundaryArg;
}


CHEAD(void) update() const {
	const BoundaryPatch& patch { this->patch() };
	Index size { patch.size() };
	Field<MeshType>& vField { field() };
	Field<MeshType>& dField { this->dField() };

	Array2Xs& vTransf { m_transformer.transform() };

	/* make sure that direction of the discharge
	 * (positive sign is inflow, negative sign is outflow)
	 * matches that of the orthogonal inner velocity,
	 * where a negative value means inflow and a positive one means outflow.
	 * This inversion is due to using outward normal vectors internally,
	 * but wanting intuitive signs for the user.
	 * std::signbit below is true when qBoundary is negative (outflow),
	 * and zero when qBoundary is positive (inflow),
	 * therefore only the case for positive qBoundary needs to be adjusted.
	 * If the velocity is zero, no error is thrown,
	 * as this is a likely initial condition. */
	if (
		(vTransf.row(0).sign() != (std::signbit( qBoundary() ) ? 1 : -1) ).any()
		&&
		!vTransf.row(0).isApproxToConstant(0, hms::epsilon)
	){
		throw BoundaryConditionError( fmt::format(
			"Diverging flow at patch \"{}\"."
			"\nSpecified discharge q={} [sqm/s] goes {}wards, "
			"but orthogonal velocity does not.",
			patch.name(), qBoundary(),
			( std::signbit( qBoundary() ) ? "out"s : "in"s )
		) );
	}

	static constexpr Index nCols { cubic::norm_noLin::bufCols };
	static_assert( decltype(m_calcBuf)::ColsAtCompileTime > nCols );
	assert( m_calcBuf.size() >= size * nCols );

	ArrayXNsMap<nCols> bufMap { m_calcBuf.data(), size, nCols };
	Array1XsMap cw { m_calcBuf.data() + bufMap.size(), 1, size };

	assert( cw.data() + cw.size() <= m_calcBuf.data() + m_calcBuf.size() );

	/* Next two instructions set up the coefficients for the cubic equation:
	 * a*x^3 + b*x^2 + c*x + d = 0,
	 * where x is the wave speed in the ghost cell cw,
	 * a = 1, 
	 * b = - 1/2 * (v_orth_inner + 2 * sqrt(g * d_inner) ),
	 * c = 0, and
	 * d = 1/2 * g * q_boundary.
	 * The original source for this is Song et al. (2011).
	 * Note especially, that q_boundary in that source is described to be
	 * always negative. Therefore, -1/2 * g * abs(q_boundary) is used here.
	 * However, there, the wave speed cw is calculated iteratively with
	 * Newton-Raphson. Here, this is avoided by using the Cardano equation for
	 * cubic polynoms.
	 * The approach here differs further from that in Song et al.,
	 * as the sign of the orthogonal velocity is removed and later added back
	 * in, to allow for a faster and real-only computation.
	 * This assumes the casus irreducibilis (pos. discriminant), in which case
	 * one solution corresponds to the subcritical case,
	 * one to the supercritical one, and one is unphysical (neg. wave speed cw).
	 * */
	/* a is 1, c is zero */
	auto b = -0.5 * (
		abs(vTransf.row(0)) +
		2*sqrt( dField.domainValues(patch)*phys::g )
	);
	auto d = -0.5 * phys::g * abs( qBoundary() );

	/* By inverting the sign of d, there is always only a single real root */
	cw = cubic::norm_noLin::root( b, d, bufMap );

	/* Buffer can now be reused, up to the address of cw */
	Array1XsMap dGhost { m_calcBuf.col(0).data(), 1, size };
	Array2XsMap vGhost { m_calcBuf.col(1).data(), 2, size };
	/* now compute d_ghost and v_ghost from cw */
	dGhost = square(cw) / phys::g;
	/* The orthogonal row is calculated, the parallel one stays the same
	 * The inverted sign is because of the convention for qBoundary, with 
	 * positive values indicating inflow and negative values indicating outflow,
	 * while outward normal vectors are used, meaning that negative values
	 * for v_ortho indicate inflow and positive ones indicate outflow. */
	vTransf.row(0) = -qBoundary() / dGhost;
	vField.ghostValues(patch) = m_transformer.retransform();

	/* all computations are required for d as well, so we set it from here */
	dField.ghostValues(patch) = dGhost;
}

CHEAD(auto) dField() const -> Field<MeshType>& {
	assert( m_dField );
	return *m_dField;
}

CHEAD(void) dField( Field<MeshType>& arg ){
	assert(arg.nComponents() == 1);
	m_dField = std::experimental::make_observer<Field<MeshType>>( &arg );
}

#undef CHEAD

#define INSTANTIATE_FOR_ALL_MESHTYPES(MeshType) \
template class OrthogonalDischarge_velocity<MeshType>;
#include "meshInstantiator.hpp"

}